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Abstract 

To the Hubbard model on a square lattice we add an interaction, W, which 
depends upon the square of a near-neighbor hopping. We use zero temper- 
ature quantum Monte Carlo simulations on lattice sizes up to 16 x 16, to 
show that at half-filling and constant value of the Hubbard repulsion, the 
interaction W triggers a quantum transition between an antiferromagnetic 
Mott insulator and a d x 2_ y 2 superconductor. With a combination of finite 
temperature quantum Monte Carlo simulations and the Maximum Entropy 
method, we study spin and charge degrees of freedom in the superconducting 
state. We give numerical evidence for the occurrence of a finite temperature 
Kosterlitz-Thouless transition to the d x i_ y i superconducting state. Above 
and below the Kosterlitz-Thouless transition temperature, Tkt, we compute 
the one-electron density of states, N(uj), the spin relaxation rate 1/Ti, as 
well as the imaginary and real part of the spin susceptibility xiQ^)- The 
spin dynamics are characterized by the vanishing of 1/Ti and divergence of 
R- e x(? = C 71 ") 71 ")) 1 ^ = 0) in the low temperature limit. As Tkt is approached 
N(uj) develops a pseudo-gap feature and below Tkt Imx(<f = (tt, 7t), u) shows 
a peak at finite frequency. 
PACS numbers: 71.27.+a, 71.30.+h, 71.10,+x 
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I. INTRODUCTION 



The motivation of this work is to examine the competition and relationship between an 
antiferromagnetic Mott insulating state and a d x 2_ y 2 superconducting state in two dimen- 
sions using numerical simulations. In particular, we are interested in the interplay of these 
two states and the remnant of the antiferromagnetic correlations in the superconducting 
state. We will see that significant antiferromagnetic fluctuations remain in the d x i_ y i su- 
perconducting state. The antiferromagnetic Mott insulator is well described by the ground 
state of the half-filled Hubbard model on a square lattice with on-site Coulomb repulsion U 
and nearest neighbor single-particle hopping t HJ. To this model, we add an extra term, W, 
which depends upon the square of the single-particle nearest-neighbor hopping. Staying at 
half-band filling and constant value of U, we have previously shown that we can generate a 
quantum transition as a function of the coupling strength, W, between an antiferromagnetic 
Mott insulating state and a d x 2_ y 2 superconducting state ||. Technically, with the system 
at half-filling, and with our particular choice of W, we encounter no fermion sign problem 
in the Quantum Monte Carlo (QMC) simulations. Large lattice sizes and low temperatures 
may thus be reached without loss of precision. While there are various ways to justify the 
form of the microscopic Hamiltonian, we view the choice of the interaction more formally as 
a means of inducing the desired quantum transition. In fact, one of the reasons for choosing 
this form is that it has a simple Hubbard-Stratonovitch representation which is useful in 
constructing the Monte Carlo simulation. Although the formalism with the W term may 
be extended to any lattice structure, we study here the two-dimensional system on a square 
lattice. 

The basic half-filled Hubbard model that we will study has the Hamiltonian 




(1) 



with the hopping kinetic energy 




(2) 
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Here, cl (c^ CT ) creates (annihilates) an electron with z-component of spin a on site i, n 7a = 
cl cq a , and 5 = ±a x ,±a y where the lattice constants. The energy will be measured 

in units of t. We consider the boundary conditions 

c I+Ld x ,a = ex P (-2?rz$/ $ ) c- )ff and c ?+L ^ j(T = c ?tT , (3) 

with $o = he/ e the flux quanta and L the linear length of the square lattice. The boundary 
conditions given by Eq. (|3|) account for a magnetic flux threading a torus on which the 
lattice is wrapped. 

The interaction that we will add has the form: 

H w = -WY,Kf (4) 

i 

with positive W. The Hamiltonian 

H = Hu + H w (5) 

has the possibility of exhibiting a quantum transition between an antiferromagnetic Mott 
insulating state and a superconducting d x 2_ y 2 phase. When W = 0, the half-filled Hubbard 
model with a finite U is known to be a Mott insulator with long-range antiferromagnetic 
order. The interaction Hw can be decomposed into the following terms: 

tt _ rr(l) , rr(2) , rr(3) „(4) 
-"VK — -"ly + -"TV + "W + n W 

h$ = -4wy d c 7a - w y cl, c UK , 

i,cr i,a,S,S' 

H$ = -W £ {^A-^,-^ +l(7 + h.c.) 

H$ = + W y (Tl ^ A1 + Tl^T^ + Tl j0 T ls>0 ) 

1,6,6* 

= ~Wy Al s Als- (6) 



Here, = cXA - , ll 1 = cl A - , , T-J = (cLcJ - + cl , cl - J /V2, and 

Aj^ = (ct^cl - — c |*j. c ]+<5t) contains single-particle terms which renormalize the 



chemical potential and allow single-particle hopping between next nearest neighbor sites. 

(2) 

H\y scatters an on-site singlet to adjacent sites. In the presence of a Hubbard interaction, 
this term should not contribute to the low energy physics since double occupancy is sup- 

(3) 

pressed. Hy/ corresponds to a triplet scattering channel. Since this term has a positive 

(3) 

coupling constant, and triplet pairing is not favored by the Hubbard interaction, Hy/ is not 
expected to be relevant for the low energy physics. Finally H$ contains the term we are 
interested in. The terms in Hyy with 5 = 5' contribute to the exchange interaction giving 

WT,(SvS U s-\ntn us ) (7) 



i.S 

r(4) 



while the terms in HyJ with 5^5' contribute to a pairing interaction. Thus in the presence 
of the on-site Hubbard repulsion, Hyf is the relevant part of Hy/ in determining the low 
energy properties. 

As previously discussed, we have chosen the form of the interaction Hyr in order of 
obtaining a model which exhibits a transition from the antiferromagnetic Mott insulating 
state to a d x 2_ y 2 superconducting state and is suitable for Monte-Carlo simulations. In 
particular integrating out the phonons from a Su-Schrieffer-Heeger |3|] term with Einstein 
oscillators: 

£ A • (Qr- Qj) (clc lr7 + h.c.) + £ (3. + q^q\ , (8 ) 

{ij),<r i \ ' 

and taking the antiadiabatic limit (M — > 0), generates H\y with W = A' D~ 1 \/2. Pairing 
mechanism along those lines were considered in [§,0. Here, however, we view this simply 
as a Hubbard Stratonovitch transformation which is useful in constructing the Monte-Carlo 
simulation. 

Our general view is that the half-filled two-dimensional Hubbard model is near various 
instabilities. In particular, Monte Carlo calculations find a divergence in the compressibility 
and an unusually large dynamical exponent at the metal insulator transition driven by the 
chemical potential. In the doped state obtained with the addition of a chemical 

potential, the leading singlet pairing interaction is found in the d x 2_ y 2 channel ll|J^]. There 



is clear evidence that the model is sensitive to the addition of a chemical potential. Here, 
we will examine the half-filled system to the interaction W. 

The organization of the text is the following. In the next section, we give a description 
of the generalizations required to implement Hw in the Projector Quantum Monte Carlo 
(PQMC) |T(]|-|r5| and finite temperature QMC algorithms P^JT5| . In the appendix, we 
discuss the construction of trial wave functions used in the PQMC. Numerical results are 
presented in Sees. (|TTl|) and (|V|). In Sec. (|T|) we concentrate on the charge degrees of 



freedom and in Sec. ([TVD on the spin degrees of freedom. The zero temperature data, from 
which we conclude the existence of a phase transition between an antiferromagnetic Mott 
insulating state and d x 2_ y 2 superconducting state has already been presented in reference 
0. At finite temperatures, we concentrate mostly on the superconducting state. We show 
the occurrence of a finite temperature Kosterlitz-Thouless transition. The one electron 
density of states, N(u), the relaxation rate l/Ti, as well as the real and imaginary part 
of the spin susceptibility are calculated in the superconducting phase above and below the 
Kosterlitz-Thouless transition temperature. In the last section, we discuss and summarize 
our results. 

II. QUANTUM MONTE CARLO ALGORITHMS 

We have applied both the PQMC and finite temperature QMC algorithms for numerical 
simulations of the Hamiltonian (j5|). The details of those algorithms in the case of the 
Hubbard model have been reviewed in the literature |2(J. In this section, we discuss the 



generalizations required for the inclusion of the Hy/ term. We will concentrate on the 
description of the PQMC algorithm. Apart from differences in the numerical stabilization of 
the algorithms, the step from the PQMC to finite temperature algorithm is straightforward 
and similar to the case of the Hubbard model. For the numerical simulations, it is convenient 
to carry out a canonical transformation, c* a = exp (2iri-^^i ■ a^j c^ a . From equation ([3]) 
it follows that the c fermionic operators satisfy periodic boundary conditions. Under this 
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canonical transformation, 

K 7 = Y (ct + ct - c 7 e l ^ ld A (9) 

where = In this section we will work in this basis, and omit the tilde on the 

fermionic operators c. The inclusion of magnetic fields in QMC methods is straightforward, 



and has been discussed in |T6 



The idea behind the PQMC algorithm is to filter out the ground state from a trial wave 
function \^t) which is required to be non-orthogonal to the ground state |\l/o) 

WW = , im (*f'° H °fP* T) do) 

The Monte Carlo evaluation of the observable O proceeds in the following way. The first 
step is to carry out a Trotter decomposition of the imaginary time propagation: 

-20H ( -AtH v -AtH ( S> -AtH ( ^ w) -ArH t 



Here, H t {Hu) denotes the kinetic (potential) term of the Hubbard model, mAr = 2B and 
Hw = E H$, H$ = -WY, K} M and [K iM , K^} = V r, r'. (12) 

n=l r=l 

The above defines the function k(n, r), n = 1 . . . n w and r = 1 . . . N/n w for an iV-site lattice. 
Hence, Hyy is given by a sum of commuting operators. The systematic error produced by 
the above Trotter decomposition is of the order At. However, provided that the operators 
Hyp, H t , H U} O, as well as the trial wave function are simultaneously real representable 
the prefactor of the systematic error proportional to At in the evaluation of the observable 
O vanishes [|18|,|l(J. As we will see, this condition is satisfied in our calculations and the 
systematic error produced by the Trotter decomposition is of order At 2 . 

Having isolated the two-body interaction terms, one may carry out Hubbard 
Stratonovitch transformations so as to express the right hand side of equation (|11|) as 
an imaginary time propagation of non-interacting electrons in an external field. For the 
Hubbard interaction, it is convenient to carry out a discrete Hubbard Stratonovitch decom- 



position |19| to obtain 



cE-p(Ei^4 (13) 

s \i,j,a / 



Here s denotes a vector of length given by the number of sites, N, of HS Ising fields and 

D U") = %cosh- 1 (Arf//2) S? a. (14) 

The constant C = exp(—ArNU/2)/2 N for the TV-site system will be dropped below. 
For the decomposition of Hffi we have used the approximate identity: 

e -A.<> =c ,,£ a (n) (f) exp ctA^ficj) + 0(Ar 4 ) (15) 
where, = (l^, l^ N /^)), lM = -2, -1, 1, 2 for r = 1 . . . N/n w and 

N/n w 



r=l 

N/n u 



The fields r\ and 7 take the values: 



7 



;±1) = 1 + VE/3, 7(±2) = 1 - v^/3 



77(±l) = ±J2(3->/6), r/(±2) = ±W2 (3 + v 7 ^) 



(17) 



The constant C" = (\/A) N l nw will also be dropped below. The matrices A^ n,r ' ) and A^ n,r '^ 
commute for all combinations of r and r'. The systematic error produced by the above 
decomposition of Hyp is of order At 4 per time slice, and hence produces an overall systematic 
error of order At 3 . Compared to the accuracy of the Trotter decomposition (of order At 2 
) this is negligible. 

The trial wave function is required to be given by a Slater determinant: 
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|*r> = |*t) ® l*r) 

N a / 

l^) = li E4,^)i° CT )- ( 18 ) 



n=l 



Here JV^ denotes the number of particles in a given spin sector, |0 CT ) is the vacuum in the 
spin-a sector and P u is an iV x Ng rectangular matrix where N is the number of sites. 

One may now integrate the fermionic degrees of freedom and define a probability distri- 
bution: 





e- 0H Oe~ 


BH 






e -2BH, 





a(f) det (M\l, s)) det (Mltf s) 

Pr(s,/) = - — — ^ =r~^- (19) 

E s -./a(/) det (MT(/, s)J det (Mi(/, s)j 

The fields s, I acquire an additional time index r : 1 . . .m, a(l) = llr,n,T7 (4™'^) an d 
M a (l, s) = P^B'ipe, 0)P a 

Bg[i e 27 e l)= f[ e D^ e A?\?)_^ e A^\n e -ArT_ (20) 

T=T1+1 

Here, 6i = nAr, 6 2 = r 2 Ar and iJ t = J2ij a ct o _T7jcj cr corresponds to the kinetic energy 
of the Hubbard model. Observables may now be evaluated by: 

J2Pr(s,l)(0) gX +0(Ar 2 ) (21) 

s,l 

For a given set of fields s, I, Wicks theorem applies and it suffices to calculate single-particle 
Green functions. They are given by: 

MAf = ( J " Vp°m°(s, ry^p^Bz^e, e))„. (22) 

Here / is the unit matrix, J" = S??. 

As mentioned previously, the Monte Carlo simulation, for the half-filled case, does not 

— * — * 

suffer from the sign problem. That Pr(s, /) is positive for all values of the fields s and I follows 
by carrying out a particle- hole transformation in say the up spin sector: cl = (—l) lx+ly hj*. 
Since the electron vacuum |0 T ) is given by |0 T ) = n^-LlO^) one has 

I, I 

/ \ N — N^ / \ 

\*b = ri \ Y,H^- i ) ix+iyp ln) ihi° m > - rf e^A io M )- (23) 

n=l \ 7 17 n=l \ • / 



The above equation defines P and Pj n denotes the complex conjugate of Pj n . Under the 
above transformation, one has: 



det (Ml f ) = e a ^ s ^det (PlB^fiG, 0)P) (24) 

where a = cosh _1 (Ar?7/2). The complex conjugation comes from the fact that if electrons 
feel a flux <3>, holes are submitted to a flux — <&. The complex conjugation changes the sign 
of the flux $. Provided that the particle number in given spin sectors satisfy N — — 
and that the trial wave function satisfies P = P^ one has: 



det (Ml -) = e a %r det ( M l ^) (25) 

— * 

from which the positivity of the probability distribution Pr(s, /) follows. We have used a 
trial wave function which satisfies the above condition for the positivity of the probability 
distribution. We furthermore require \^t) to be a spin singlet: 

E-Mj|*t> = (26) 

— * — » 

where is the spin operator on site %. An explicit construction of trial wave functions 
showing no sign problem, being spin singlet and if necessary real representable in real space 
is given in the appendix. 

The sum over the Hubbard Stratonovitch fields is carried out with Monte Carlo methods. 
One sweeps sequentially through space time lattice and proposes single spin flip updates. 
The ratio of new to old probabilities under a local change /[™' r - ) — ■> Z^ n,r ) is best calculated in a 
basis where the matrix J^f^ is diagonal. Since in this basis A^-f^ has only two non-vanishing 
eigenvalues, the ratio of probabilities involves the knowledge of four Green functions per spin 
sector. Note that under a proposed local move s^ T — > the ratio of probabilities is given by 
a single Green function per spin sector. In case of acceptance of the spin flip, the required 
updating is two times more expensive in the case of an Z[ n,r ) — > /^ n,r ) move than for an 
s^ T — > s'^ r move. The numerical stabilization of the algorithm is identical to that used for 
Hubbard calculations. This concludes our description of the PQMC algorithm. 
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The finite temperature QMC algorithm, is a grand canonical simulation which evaluates: 

Tr (e- pH d) 

<°> = ^(Fmr < 27 > 

where the trace runs over the Fock space and f3 denotes the inverse temperature. The finite 
temperature QMC algorithm is conceptually similar to that of the PQMC |[20|| . Both the 
PQMC [17] and finite temperature QMC algorithms may be used to calculate imaginary 
time displaced correlation functions. The dynamical results presented here stem from the 
use of the Classic Maximum Entropy method |2l|j2^| to analytically continue imaginary time 
data produced by the finite temperature QMC algorithm. 

We are now in a position to test and compare the PQMC and finite temperature algo- 
rithms. As mentioned previously, the systematic error produced by the Trotter decomposi- 
tion is of the order At 2 provided that the operators , H t , Hjj, O, as well as the trial 
wave function \^t) are simultaneously real represent able. This may be checked explicitly 
by calculating the energy with the PQMC algorithm at $/$ = 0.25, Z7/t = 4, W/t = 0.35 
on a 4 x 4 lattice with 20t = 1 for various values of At. The results are plotted in Fig. 
[1]. They follow well an a + bAr 2 + cAr 3 form. Note that the At 3 term, originates from 
the approximate Hubbard Stratonovitch decomposition of H\y. We have carried out most 
of our PQMC simulations at At = 0.0625t which is sufficiently small for our purposes. 

Fig. shows plots of the total energy E as well as the spin structure factor at Q = (n, 7r), 
S(Q) = | J2r elQr (SrStf) as obtained with the PQMC and finite temperature algorithms. For 
the finite temperature data, we set (3 = 20 and measure the observables with equation (^7|) . 
The same observables are evaluated with equation fl2~l|) and the above described trial wave 
function. We set At = 0.1 and in the case of the PQMC algorithm measure observables on 
the ten central time slices. For observables which do not commute with the Hamiltonian, 
this yields an effective projection parameter O e ff = @ — 0.51 We consider a 6 x 6 lattice 
at U/t — 4, W/t = 0.35 and $ = 0. The trial wave function used here, is constructed as 
shown in Eq. ( |A4[ ) and the orthogonal transformation used is obtained by diagonalizing 
the Hamiltonian of equation (|A6|) at a numerically infinitesimal value of S. We have used 
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this trial wave function for all our calculations at $ = 0. In the limit of large 0t, both 
the finite temperature and PQMC results converge to the same value within the error bars. 
We observe that ground state properties are obtained more efficiently within the PQMC 
approach. 



III. CHARGE DEGREES OF FREEDOM 

In this section, we consider the charge degrees of freedom at zero and finite temperatures. 
We start by showing that at zero temperature, half-band filling, and constant value of the 
Hubbard repulsion, U/t = 4, W/t triggers a transition between an insulator and a d x 2_ y 2 
superconductor. The transition is found to occur, at W c /t ~ 0.3. We then consider the 
superconducting state at W/t = 0.35 and give numerical evidence for the occurrence of 
a finite temperature Kosterlitz-Thouless transition. The temperature dependence of the 
one-electron density of states, N(u), is analyzed. 



A. Zero temperature 

An efficient and general method to distinguish between an insulator and superconductor 
is to compute the ground state energy as a function of a twist $ in the boundary condition 
along one lattice direction: E (Q). In the case of an insulator, the wave function is localized 
and hence, an exponential decay of AEq(&) = Eo(Q) — Eq(^q/2) as a function of lattice size 
is expected ||23|| . In the spin density wave (SDW) approximation for the half-filled Hubbard 



model, one obtains Ai?o($) = exp (— L/^) where £ is the localization length of the 

wavefunction. On the other hand, for a superconductor (i.e off-diagonal long-range order 
(ODLRO)), AEq(<&) shows anomalous flux quantization: AE (&) is a periodic function of $ 
with period $o/2 and a non vanishing energy barrier is to be found between the flux minima 
p4|- p6| , p!6| so that AEq(^q/A) remains finite as L —>■ oo. In general, this functional form of 
AEq(&) should occur only in the thermodynamic limit. 
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Fig. | shows A£ ($) at W/t = 0.1 and W/t = 0.5. At W/t = 0.1, and as expected for 
an insulator, AE (&) is suppressed with growing lattice sizes for all values of $. In contrast, 
at W/t = 0.5 one notes that for values of 1/4 < < 1/2, A£ , ($o) is very stable against 
growing lattice sizes. On the other hand, at $ = 0, AEq(&) decreases with growing lattice 
sizes. We interpret the data as a finite size indication of anomalous flux quantization. At 



the end of Sec. ( |IV A| ) we will argue that this slow size scaling of AE (&) in the vicinity 
of $ = is related to the nature of the spin degrees of freedom in the superconducting 
state. Very similar results were obtained in the case of the repulsive and attractive Hubbard 
models at half-filling |26|,|]1|. 



To best distinguish between the insulating and superconducting states, we consider the 
size scaling of the quantity A£ , ( < ^ ) o/4) for various values of W/t as in Fig. f|a. One observes 
a change in the size-scaling of A_E' ($ /4) as W/t decreases from W/t = 0.5 to W/t = 0.22. 
From these measurements, we estimate that the change occurs in the vicinity of W/t = 0.3. 
For values of W/t < 0.3, A£ , ($ /4) is consistent with the SDW form whereas for W/t > 0.33 
A£'o($o/4) may be extrapolated to a finite value with finite size corrections proportional 
to 1/L. The extrapolated value of Ai? ($o/4) versus W/t is plotted in Fig. f|b and the 
quantum transition between a Mott insulator and superconductor occurs at W/t ~ 0.3. 

In order to determine the symmetry of the order parameter in the superconducting state, 
we have calculated pair- field correlations in the s and d x 2_ y 2 channels: 

PdM = (Ai )S (r)A d , s (0)) (28) 

with 

Ai, s (rO = £ UM°4A + l-« ( 29 ) 

Here, f s (5) = 1 and f d (6) = 1(-1) for 6 = ±a x (±a y ). Fig. § shows plots of P d>s {L/2, L/2) 
for W/t = 0.6 where the system is superconducting and W/t = 0.1 where it is not. At 
W/t = 0.6, the dominant signal at long distances (L = 16) is obtained in the d x 2_ y 2 channel. 
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At the mean-field level, the symmetry of the order parameter will determine the func- 
tional form of the single-particle occupation number, n(k). For a d x 2_ y 2 superconductor the 
BCS result yields: 



n(k) = 1 , Cfc (30) 

lAl + e? 



k k 

where corresponds to the 

single-particle dispersion relation and = Aa (cos(fc x ) — cos(k y )). As is apparent from 
Eq. (|30|), in the = fc(l, 1) direction where the d x i_ y i gap vanishes n(k) shows a jump 
at the Fermi energy, whereas in the k = k(l,0) direction n{k) is a smooth function of k. 
Precisely this behavior in n(k) may be detected in the QMC data at W/t = 0.6 as shown in 
Fig||a. For comparison, we have plotted n(k) at W = where it is expected to scale to a 
smooth function in the thermodynamic limit (see Fig||b). 

B. Finite temperature: W/t = 0.35 

For values of W/t > W c /t ~ 0.3, one expects the occurrence of a finite temperature 
Kosterlitz-Thouless transition. To detect this transition we consider the quantity: 

AE ($, T) = E ($ = $ /4, T)-E($ = $ /2, T) . (31) 



In the framework of a Kosterlitz-Thouless p7| , p8| transition one expects in the thermody- 
namic limit 

AE ($, T) ~ 8{T — T KT ) (32) 

where Tkt is the Kosterlitz-Thouless transition temperature. The above follows by consid- 
ering the Free energy = — jjj In Tre _/3H ^*' ) . Expanding F around $ = $o/2 yields: 

= F&/2) + i f^)^ + O ((^^) 4 I where 
£> s (/?) = - (y) 2 (< K x >+J^dr< J x (r)J x (0) > 
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i 

K X = ^-2W (jf • j, ,. f + W {K r +K r+Sx ,jf} 

i 

j? = — i ( ct c?,~ — ct _ Ct> ) and 

J i ^ \ i,a i+a.x,cr i+a x ,cr ) 

k? = - V (ct c*, n „ + ct . c 7 ) . (33) 

The boundary conditions for the calculation of the superfluid density D s satisfies: Cf +LSx a = 
—c?„ and c?, r , „ = c*„. Hence, 

Z,f7 l-\-L/CLy,(7 2,(7 ' 

A£ ($, T) = A (0 F (*) - /3F($ /2)) ~ 1 ^~^° /2 ) - T A Ds ") (34 ) 



At TicT) -Ds is expected to show a universal jump |28f so that -Jf-Ds behaves like a Dirac 



^-function and equation (|32|) follows. Such signatures of the Kosterlitz-Thouless transition 
have already been observed in classical 2J| and quantum XY models [|30| as well as in the 



attractive Hubbard model |31|Jll 

Figure [7|a plots AE(<&,T) at W/t = 0.35 and U/t — 4 for lattice sizes ranging from 
L = 6 to L = 10. For each considered lattice size, AE(<b,T) shows a maximum at finite 
temperature. For a given lattice size, we define Tkt by this maximum and obtain: Tkt ~ 
0.95, 0.75, 0.65 for the lattice sizes L = 6, 8, 10 respectively. We note that the scaling of 
the peak height is influenced by the size dependence of Tkt due to the linear T term 
in front of -§fD s (see equation (B3I)). The extrapolation of Tkt to the thermodynamic 



dT J 

limit is hard. However, we recall that the size scaling of the zero temperature energy 
difference AE (§ /4) followed well a + b/L form. Since AE (^ /A) is proportional to the 
zero temperature superfluid density, it is plausible to assume that the size dependence of 
Tkt equally follows an a + b/L form. The above three finite size values of Tkt are consistent 
with this Ansatz, and we obtain 

T KT ~ 0.2* + 4.5*/L at W/t = 0.35. (35) 



Figure |7|b plots the vertex contribution to the equal time pair-field correlations [32| both 
in the s and d x 2_ y 2 channels for the L = 6 and L = 8 lattices at the largest distance: 
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— * 

R = (L/2, L/2). This quantity is given by: 

PMfi = p*.W - E hs&hsfi) ((4 lCT ^ iCT )(c!. +5>a c QV(J ) + (ci a c^)(^ +l _ a c s , t _ a )) . 

a,5,P 

(36) 

Per definition, P^ s (r) = in the absence of interactions. An enhancement of the <i-wave 
signal at an energy scale set by Tkt is observed. In contrast, the s-wave response vanishes 
for all considered temperatures. Contrary to the AE(<&,T) data, those calculations were 
carried out at $ = 0. 

We now consider the one-electron density of states: N(u>) = j^J2kA(k,u)) with 

G{k ^ = -J-J"TT^ A{k ^ (37) 

where G(k, r) = J2 a ( c k o-( r ) c s )' r > and cg CT (r) = e TH c^ cr e~ TH . We have used the Classic 
Maximum Entropy method to obtain dynamical data from imaginary time Green functions. 
We chose a flat default model for N(u) and enforced the following sum rules: 



duN(uj) = 2vr 

J duN(u;)ujtanh({3uj/2) = ^£<[4 )(T , [#,%J]> (38) 

The calculations were carried out at $ = $o/2. The data at temperatures lesser and greater 
than Tkt is plotted in Fig. ||. We consider 6 x 6 to 10 x 10 lattices. For all three considered 
lattices sizes, one sees the onset of a pseudo-gap in the temperature range of T KT . At lower 
temperatures, than shown in Fig. || N(u) suffers form size effects. The onset of finite size 
effects on small lattices (i.e. L = 6 and L = 8 ) coincides approximately with the magnetic 
scale J ~ 0.5t, which will be determined in the next section. As the lattice size increases, one 
can go below the J scale without noticing an anomaly in N(u). This may be seen explicitly 
in Fig. |8](m) where T = 0.33t < J. We may estimate the value of the superconducting 
gap A sc by the peak position in N(u), at the lowest temperature scale presented in Fig. 
^|. For the three considered lattice sizes, the data is consistent with A sc /Tkt ~ 2.5. The 
peak position decreases as a function of decreasing temperature. On the 10 x 10 lattice, we 
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have looked more systematically at temperatures above Tkt- The data shows a pseudo-gap 
feature above Tkt- It is not clear if this aspect of the data will survive in the thermodynamic 
limit. 

IV. SPIN DEGREES OF FREEDOM 

As in the previous section, we first concentrate on the zero temperature spin correlations 
and show that the long-range antiferromagnetic order disappears at W c /t ~ 0.3. We then 
study finite temperature spin dynamics in the superconducting phase at W/t = 0.35. 

A. Zero temperature 

To detect the existence of long-range magnetic order, we compute the real space spin-spin 
correlations: 

S(r) = ^(S(r)S(0)) (39) 

— * 

where S(r) denotes the spin operator on site r. For values of W/t < 0.3 and lattice sizes 
ranging from L = 4 to L = 12, S(L/2, L/2), may be fitted to a 1/L form and scales to a finite 
value, as shown in Fig. ||a We therefore conclude that long-range antiferromagnetic order is 
present for W/t < 0.3. The associated staggered moment, m = lim£_ +00 y / 3S'(L/2, L/2), is 
plotted in Fig. ^o. The data is consistent with a continuous decay of m as W/t increases 
towards 0.3. At W/t = 0.3, we were unable to distinguish m from zero within our statistical 
uncertainty. Hence, we conclude that long-range antiferromagnetic order vanishes at W/t ~ 
0.3. Therefore, within our numerical resolution, the antiferromagnetic transition point is 
not separated from the superconductor-insulator transition point. Well within the d x 2_ y 2 
superconducting phase the spin-spin correlations remain sizable. In fact, at W/t = 0.6, 
lattice sizes ranging from L = 4 to L = 16, S(L/2, L/2) scales as L~ a with a — 1.16 ± 0.01 
as shown in the inset of Fig. |9|a. This power-law decay of the spin-spin correlations in 
the superconducting state arise because of the nodes in the d x 2_ y 2 gap. Following Ref. 
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|33| , one can approximate the spin susceptibility, x(<?> *^m), in the superconducting state 
by inserting the irreducible BCS spin susceptibility, \o (<fj ^f m ), in the random phase 
approximation (RPA) form of x(q, iv m ): XRPA(q,iuj m ) = Xo CS (& iu m )/(l-Uxo CS (<l, iu m )). 
Here, u; m corresponds to Matsubara frequencies. Within this approximation and at half- 
band filling, the spin-spin correlations for a d x 2_ y 2 superconducting order parameter show a 
power- law decay: S RPA (L/2, L/2) ~ L~ 3 - 5 . In contrast, the spin-spin correlations in an s- 
wave superconductor decay exponentially. Thus, the fact that spin-spin correlations remain 
critical in the superconducting state, may be attributed to the node and phase factors of 
the d x 2_ y 2 superconducting gap. 

The very slow decay of the spin-spin correlations in the superconducting state revealed by 
the QMC data shows the extreme compatibility between a d x 2_ y 2 superconductor and sub- 
stantial short-range antiferromagnetic spin-spin correlations. In the superconducting state, 
the spin stiffness vanishes since spin-spin correlations decay more rapidly than cos(Qf)/\r\, 
where Q = (tt,tv). However, since S(f) decays slower than cos(Qr)/\f\ 2 and the dimension 
is equal to two, S(Q) ~ ^^e 1 ^ S{f) diverges. 

At this point, we can offer an explanation of why it is convenient to study charge degrees 
of freedom at <3>/$o = 0.5 and spin degrees of freedom at $ = 0. We have just mentioned that 
in the superconducting state, the spin-spin correlations decay as l/|r| a where a is slightly 
larger than one ( a. ~ 1.16 at W/t = 0.6). On a finite lattice, it is hard to distinguish between 
true long-range order and this slow decay. Choosing anti-periodic boundary conditions in 
one lattice direction and periodic boundary condition in the other (i.e. $/$o = 0.5 ), renders 
the ground state of the half-filled non-interacting system (U = W = 0) non-degenerate. In 
a weak coupling approach, one may understand that this choice of boundary conditions 
introduces frustration in antiferromagnetic correlations since on a finite lattice XoiQ,^ — 0) 
takes a finite value in the zero temperature limit. In contrast, and on a finite lattice with 
periodic boundary conditions (i.e. $ = 0) xoiQ,^ — 0) diverges in the zero temperature 
limit. Here, xo denotes the spin susceptibility at U = W = 0. This frustration of the spin 
degrees of freedom introduced by the boundary condition, minimizes the size dependence 
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of charge degrees of freedom. This pathology may be seen explicitly in Fig. where the 
superfluid density should in principle be given by the curvature of AE (Q) at $ = or 
$ = $o/2. At $ = 0, the size scaling of this quantity is hard to extract from simulations up 
to L = 12. 



B. Finite temperature: W/t = 0.35 



We first consider the static spin susceptibility, 

ri 3 



Rex(<?> = 0) = / (m z (q,T)m z (-q,0)) 
Jo 



m z (q, t) = e rH m z (q) e~ rH and m z = (n r - T - n? A ) . (40) 



Before considering the superconducting phase, we briefly summarize the temperature depen- 
dence of Rex (q, oj = 0) at wave vectors q = and q = Q = (tt, tt) for the Hubbard model at 



U/t = 4. The results are plotted in Fig. IT^. One expects the physics of the half- filled Hub- 
bard model to be well described by a Heisenberg model in the renormalized classical regime 
|34|| . Hence, an exponential divergence of R,ex{Q, = 0) as a function of temperature is ex- 



pected. The data (see Fig. [I~0|b) shows a sharp increase at a temperature scale T ~ 0.25. We 
use this criterion to define the magnetic energy scale J ~ 0.25 at W/t = 0. At temperatures 
T > J the uniform static susceptibility Rex (q — 0, uj — 0) increases monotonously with de- 
creasing temperature. In this temperature regime, the overall scale of Rex {q = 0,u = 0), 
is greater than its free electron value. This enhancement can be understood by the Stoner 
enhancement in the RPA approximation. In the low temperature limit Rex (q — 0, o> = 0) 
scales to a finite value due to the existence of gapless spin-wave excitations. One notices 
that Rex (<f = 0, u = 0) shows a maximum at approximately T ~ J. 

At W/t = 0.35 Rex(Q, oj = 0) diverges in the low temperature limit (see Fig. |ll|b). This 
divergence is weaker than at W/t = 0, and reflects the slow decay of the real space spin- 
spin correlations. In the superconducting state the spin stiffness vanishes, and one expects a 
power-law divergence of Rex(Q, oj = 0) as a function of temperature. From this data, we may 
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identify a magnetic energy scale J ~ 0.5£ at W/t = 0.35. Being a local quantity, J is to first 
approximation insensitive to lattice size. That J is greater at W/t = 0.35 than at W/t = 
may be seen by inspecting the nature of the added interaction H\y- In the introduction, we 
have argued that at finite values of U/t, the relevant physics contained in the H\y term is 
given by H$ in equation (^). As noted, H$ contain terms of the form — WA| Aj 5 which 
may be written as W (jf? ■ S^ + g — \n^ + ^j . This term explicitly enhances the value of J. 
At W/t = 0.35, and over the entire considered temperature region, Rex (q = 0,u = 0) is 
suppressed in comparison to its free electron value (see Fig. |ll|a). We interpret this overall 
suppression as a signature of pairing fluctuations. For our two largest lattice sizes, L = 6 
and L = 8, Rex (<f = 0, u; = 0) shows a maximum at T ~ 0.5t which coincides with our 
previously determined J scale. The fact that this maximum does not scale with system size 
(at least for our two largest lattices) makes us believe that it is related to the J scale and 
not to Tkt- In the low temperature limit, and as expected for a superconducting state, 
the QMC data is consistent with the vanishing of Rex (<f = 0, to = 0). We note that in BCS 
theory and for a d x 2_ y 2 order parameter, one obtains: Rex (<f = 0, uo = 0) ~ T in the low 
temperature limit. 

T KT is determined by using different boundary conditions (see equation (|3~TD) than 
Rex(q, u = 0), which we calculate at $ = (see equation @). Hence care has to be 
taken when comparing those two quantities. Strictly speaking, one should extrapolate the 
two data sets to the thermodynamic limit and only then compare. This is especially true in 
our case, since the size effects are strong in the determination of T KT (see equation fl35D). 
We have seen that Tkt scales to 0.2t in the thermodynamic limit. At this temperature scale, 



no evident anomaly appears with growing lattice size in the QMC data in Fig. |TT 
We now consider the dynamical structure factor: 

S(q,u)= 1 _ e _ Pw (41) 

where 

roo 

X (q, u) = -i dte^+^(K (q, t) , m z (-q, 0)]) and m z (q, t) = e itH m z (q) e~ itH . (42) 
Jo 
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In the above, we have set H = 1 and 5 is a positive infinitesimal number. lmx{q, oS) was 
obtained by analytically continuing the imaginary time QMC data with the use of the Classic 
Maximum Entropy method. We have used a flat default model for S(q,u) and imposed the 
sum rules: 

lmx(q, u) 



J-oo UJ 

dulmx(q, co>)coth (flu/2) = 2n(m z (q) m z (—q)) and 
dulmx(q,uj)uj = n([m z (-q) , [H,m z (q)\\). (43) 



oo 
oo 



-oo 
oo 



Fig. [12] plots S(Q, u) as a function of temperature at W/t = 0.35 on 6 x 6 and 8x8 lattices. 
For comparison, we have plotted at the lowest considered temperature, (3 = O.lt, S(Q,u) 
for the Hubbard model (i.e. W/t = 0). S(q,ui) satisfies the sum rule: 

/oo 
duS(q,u) = Tc(m z (q) m z (-q)) = nS(q). (44) 
-oo 

At W/t = 0, S(Q) diverges as L 2 in the zero temperature limit. This weight is centered 
at uj = 0. At W/t = 0.35, and at a scale set by J ~ 0.5t, we observe a buildup of 
spectral weight centered around u — 0. At a lower temperature scale, which we will identify 
below, spectral weight is shifted from lower frequencies to form a peak at finite frequency. 
Defining u by the peak value of S(Q, uj) at the lowest considered temperature, we see that 
ujq/Txt ~ 0.4. This relation is valid for both considered lattice sizes. Since the equal time 
spin-spin correlations in the superconducting state decay slower than cos(Qr)/\r\ 2 at zero 
temperature, the spectral weight under the peak at ujq diverges in the thermodynamic and 
low temperature limits. 

The relaxation rate we consider is defined by, 

1 -Js^E 52 ^. («) 



1 q 

where a ^-independent nuclear form factor is assumed. Fig. |13| plots \/T\ both in the 
superconducting and antiferromagnetic insulating state. The error bars are obtained in the 
following way. For given imaginary time data and covariance matrix, we transform the data 
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into a basis where the covariance matrix is diagonal. In this basis, we add independent 
gaussian noise to each data point. The width of each gaussian distribution is determined 
by the diagonal of the covariance matrix. We carry out the analytical continuation and 



repeat the procedure for several realization of the gaussian noise. The error bars in Fig. |13 
correspond to the variance of the so obtained values of 1/T\. 1/Ti is a delicate quantity to 
compute with QMC methods. To check the reliability of our calculations we first consider 
the W/t = case. At temperatures T < J as determined from R,ex{Q,uJ = 0), an increase 
in 1/Ti is observed. This is an expected feature since in the renormalized classical regime 



of quantum antiferromagnets, l/Ti diverges exponentially with decreasing temperature |34 
In the superconducting state at W/t = 0.35, we plot l/Ti as a function of T/T KT (see Fig. 
P~3|a) . The overall scale of l/Ti is reduced in the superconducting state as compared to 
the antiferromagnetic Mott insulating state. For both considered lattice sizes, l/Ti shows 
a maximum at T ~ 0.257Vt which allows us to define a cross-over temperature scale: 
T^/Tkt ~ 0.25. At temperatures lower than T£ r , the QMC data is consistent with the 
vanishing of the relaxation rate l/Ti. In BCS theory with a d x 2_ y 2 order parameter, l/Ti ~ 
T 3 in the low temperature limit. At temperatures above Tf , l/Ti grows with decreasing 
temperature. One expects l/Ti to start increasing at the magnetic scale J as determined 
from Rex(Q,u = 0). Since in Fig. |l3| a we set the temperature scale by Tkt large size 
effects are now present in J/T^t which takes the values 0.5 (0.65) for the L = 6 (L = 8) 
lattice and scales to J /Tkt ~ 2.5 in the thermodynamic limit. The so determined cross-over 
temperature coincides well with the temperature scale at which spectral weight in S(Q,u) 
is shifted from low frequencies to form the peak centered at ujq (see Fig. [12]). 

V. SUMMARY AND CONCLUSIONS 



As schematically illustrated in Fig. |H|a our zero temperature results at U/t = 4 are best 



summarized in the W — fi plane where fi denotes the chemical potential. We have carried 
out our simulations at half-filling, /i = 0, and shown the occurrence of a quantum transition 
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between an antiferromagnetic Mott insulator and d x 2_ y 2 superconductor as a function of 
W/t. Within our numerical accuracy, the quantum transition occurs simultaneously in both 
the charge and spin degrees of freedom at W c /t ~ 0.3. Although the spin stiffness vanishes in 
the superconducting state, the antiferromagnetic spin correlations decay slower than l/|r| 2 
where \r\ denotes the distance. In two dimensions, this slow decay leads to a divergence 
of the equal-time spin structure factor at wave vector Q = (71,71). This feature shows the 
extreme compatibility of antiferromagnetic fluctuations and d x 2_ y 2 superconductivity. To a 
first approximation, the exponent of the power-law decay of the spin-spin correlations varies 
with the coupling strength W/t. At W/t = 0.6 the QMC data for lattice sizes ranging 
from L = 4toL = 16is consistent with real space spin-spin correlations of the form 
cos(Qr)/|r*|~ a with a ~ 1.16. The QMC data suggests that a decreases continuously to 
a = 1 as the coupling W/t approaches W c /t. Hence, W c /t ~ 0.3 terminates a critical line 
in the schematic phase diagram drawn in Fig. Ok. 



As a natural consequence of Fig. |TJ]a we expect at W < W c a doping induced quan- 
tum transition between an antiferromagnetic Mott insulator and d x 2_ y 2 superconductor. As 
mentioned in the introduction, the anomalous compressibility and unusually large dynami- 
cal exponent observed numerically at the filling controlled metal-insulator transition |6]-§. 
enhances the sensitivity of the system to two-particle processes as generically contained in 
the form H^. Work along those lines, for the filling controlled transition, is under progress. 

To determine the order of the phase transition at W c /t is presently beyond the reach 
of our numerical calculations. At the mean field level, the transition is expected to be of 
first order since the phases have different broken symmetries. In this scenario, the staggered 
moment would be expected to show a jump at the critical coupling constant. Within our 
numerical accuracy, we do not observe such a feature, and the possibility of a continuous 
phase transition remains open. An example of a continuous phase transition between two 
broken symmetry states on a square lattice is found in the case of quantum antiferromag- 
nets where the Berry phase is a dangerously irrelevant operator which leads to spin-Peierls 



ordering (i.e. broken lattice symmetry) in the disordered state for half-integer spin |35H37 

22 



In the k — type BEDT — TTF compounds, the antiferromagnetic and superconducting 
phases are adjacent to each other in the plane of temperature and either pressure, anion sub- 
stitution or deuteration of hydrogen atoms |38|-f40|. Thus both changes in the bandwidth 



and in the interaction strength have been studied. The insulating phase of the n— type 
BEDT — TTF compounds is a correlated insulator in the sense that band-structure calcu- 
lations predict a metal JO] . A dimer model has been proposed to account for the magnetic 



insulating phase ||38|| . Here, a pair of BEDT — TTF molecules carries a charge of unity, and 
constitute a single site in terms of a Hubbard model. The on-site Hubbard interaction Udimer 
depends upon the intra-dimer hopping, tdimer- In the limit of large Coulomb repulsion per 
BEDT — TTF molecule, Maimer ~ ^dimer- Because of the layered structure of this compound, 
our two-dimensional model may offer a simplified description of the system. In those com- 
pounds, the direct transition line between the antiferromagnetic and superconducting phases 
appears to extend to finite temperatures, thus implying a first order phase transition. A de- 
tailed comparison between theory and experiment is beyond the scope of this work. However, 
an interesting point, is that in the superconducting phase both our results and experimental 
results show a common feature: the peak temperature in Rex(? = 0, u = 0) is higher than 
the crossover temperature in 1/Ti (i.e. Tf 7 * in our notation). It would be interesting to see 
whether the antiferromagnetic fluctuations are robust (i.e. divergence of R,ex{Q,uJ = 0)) 
in the superconducting phase. We also note that the 13 C nuclear spin relaxation rate in 
the superconducting phase of the n— type BEDT — TTF compounds has been reported to 
follow a T 3 law ||2| which would be consistent with a d x 2_ y 2 order parameter. 

We have studied spin and charge degrees of freedom at finite temperatures in the super- 



conducting phase, at coupling strength W/t = 0.35. As schematically drawn in Fig. |14| b we 
expect the occurrence of a finite temperature Kosterlitz-Thouless transition. We give numer- 
ical evidence that such a transition indeed occurs. After extrapolation to the thermodynamic 
limit we estimate Tkt ~ 0.2i at W/t = 0.35. In the vicinity of the Kosterlitz-Thouless tran- 
sition temperature a pseudo-gap in the one-electron density of states, N(u), appears. On 
our finite sized lattices, the pseudo-gap feature appears above Tkt- For the three considered 
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lattice, L = 6, L = 8 and L = 10, the estimated superconducting gap scales as T KT and 
satisfies A SC /T KT ~ 2.5. 

In the superconducting phase, at W/t = 0.35, Rex(Q,u = 0) diverges in the low tem- 
perature and thermodynamic limits. This divergence as a function of temperature is slower 
than in the antiferromagnetic Mott insulating state at W/t = 0, and is expected to follow a 
power-law. We may estimate a magnetic scale by the temperature at which Rex(Q, oj — 0) 
starts to increase. At W/t = 0.35, we obtain J ~ 0.5t The QMC data is consistent with 
the vanishing of Rex(<f = 0, oj = 0) in the zero temperature and thermodynamic limits. 
R e x(<? — 0, oj — 0) shows a maximum at T ~ J. A similar although stronger feature is 
seen in the Hubbard model (i.e. W/t = 0). From the relaxation rate 1/Ti, we estimate a 
cross-over temperature, Tf r below which 1/Ti decreases. This temperature scale satisfies 
T± r /T KT ~ 0.25. In the temperature range < T < J, 1/Ti grows with decreasing tem- 
perature. At the temperature scale set by T{ r spectral weight in S(Q, oj) is shifted from low 
frequencies to form a peak at oj . The location of the peak, oj , scales with T KT as a function 
of system size and satisfies oj /T kt ~ 0.4. The spectral weight under this peak diverges in 
the zero temperature and thermodynamic limits. 

The superconducting state at W/t = 0.35 is thus best characterized by i) the divergence 
of Rex(Q,uj = 0), but the vanishing of 1/Ti in the low temperature and thermodynamic 
limits and ii) oj /A sc ~ 0.15 << 1. This behavior is non-trivial and points out to the 
extreme compatibility of d x 2_ y 2 superconductivity and antiferromagnetic spin fluctuations. 

To conclude, the here presented model, shows a rich phase diagram. The doping of the 
model offers the possibility of reproducing some of the features of the cuprates at an energy 
scale accessible to numerical simulations. 

* Future address: Institut fur Theoretische Physik III, Universitat Stuttgart, Pfaffenwaldring 
57, D-70550 Stuttgart, Germany. 
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APPENDIX A: TRIAL WAVE FUNCTIONS FOR THE PQMC 

In this appendix we show how to construct trial wave functions which generate no sign 
problem if + iV^ = N, and in the special case, = = N/2 are spin singlets. We 
provide explicit forms of the trial wave functions used in this work. Our starting point is 
the non-interaction Hamiltonian: 

#o = E4% c j- ( A1 ) 

hi 

Ho corresponds to the Hamiltonian (|5|) in the limit U = W = 0. We have omitted the spin 
index, and consider the boundary condition given by (H). Let U be the unitary transfor- 
mation which diagonalizes T: U^TU = diag (Ai . . . A^v) . We define the hole operators by 
ht = (— lVct as well as 

1% = E Vlfi ^d r, £ = £ (A2) 

i i 

We have used the notation (—1)' = (— V) lx+l y . The above operators satisfy 
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[#°, 7S = A^l and [ff 0l 4] = -A^t. (A3) 
We set the trial wave function in the up spin sector to be the ground state of H : 
\^t) — HkeA'kl^)' wnere |0) is ^ ne electron vacuum and A is set of k points defined by 
the requirement that \^t) is the ground state of Hq. By construction \^ T ) = II^ ^lO' 1 ), 
where \0 h ) is the hole vacuum and A is the complement of A. We may hence define: 

Ph =U VL ^l---^ keA 

Pj =(-lfU u i:l...N, keA. (A4) 



With this choice of P a no sign problem occurs (see Eqs.(|T8|), (|23|), (|24| ) and (p5|)). We note 
that the only restriction on the particle number is A^ + A^ = A. 

In the special case A^ = = A/2 the trial wave function \^t) — I^t) ® \^t) satisfies: 

£<*t$- SjI^t) = 1 -J2P(k)P(k) + (P(k) - l) (P(k) - l) where 

hi k 

P{k) = (^1^1^) and P(k) = (VtilteH)- (A5) 

If the ground state of H with A/2 particles per spin sector is non-degenerate, then 
J2i j{^t\S^ ■ Sj\^t) vanishes, independently of the choice of the unitary transformation U. 
For non-zero values of the flux $ (i.e. $ 7^ 2nn^o, n being an integer), this condition is 
satisfied, and hence a trial wave function showing no sign problem and satisfying (|26|), may 
easily be constructed numerically. 

At $ = 0, the ground state of H with A/2 particles per spin sector is degenerate, and 
the total spin of the trial wave function ill depend on the choice of the unitary transformation 
U . In principle, one can use the trial wave function obtained at an infinitesimal value of $ 
to lift the degenerancy. This procedure would force us to work with complex numbers which 
enhances the CPU time by a factor three to four. To easily circumvent those problems we 
define a new Hamiltonian 

Ho = -f (1 - *) £ (c^ + h.c.) -t(l-5)^ £ + h.c.) 

i i=(2i x +l,i y ) 

-Ki + s)^ E (4 c ^ + H ( A6 ) 
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with periodic boundary conditions in both lattice directions. For infinitesimal values of 5, 
the ground state of H with N/2 particles is non-degenerate. Diagonalizing H generates an 
orthogonal matrix U from which one may construct the trial wave function. Both for the 
here presented model (see Fig. |^) and for the half-filled Hubbard model |T7[], the above trial 
wave functions provide quick convergence to the ground state. The reason lies in the fact 
that they are spin singlets and hence orthogonal to spin-wave excitations which constitute 
the low energy part of the spectrum at least in the case of the half-filled Hubbard model. 
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FIGURES 




FIG. 1. Total energy as obtained with the PQMC at Qt = 1 as a function of At. The solid 
line is a least square fit to the form a + bAr 2 + cAr 3 with a = -81.146 ± 0.015, b = -153 ± 3 and 
c = 334 ±15. The energy unit is set by t = 1. 
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6x6, (n) = 1, U/t = 4, W/t = 0.35, $ = 




FIG. 2. (a) Total energy E and (b) spin structure factor at Q = (ir,ir) as obtained from the 
finite temperature QMC algorithm at f3 = 29 and from the PQMC algorithm. 
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A£ ($) U/t = 4, (n) = 1 
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FIG. 3. AE {<$>) 



= E (<Z>) - E {<S> /2) at (a) W/t 
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= 0.1 (b) and W/t = 0.5 
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FIG. 4. (a) AE (<S> /4) = E (<£> /4) - E (<S> /2) versus l/L for several values of W/t. For 
W/t < 0.3 the solid lines correspond to a least square fit of the data to the SDW form: L exp(— L/£). 
For W/t > 0.3 the QMC data is compatible with a l/L scaling to a finite constant. The solid lines 
are least square fit to this form, (b) Extrapolated value of A£'o(^ > o/4) versus W/t. The solid line 
is a guide to the eye. 
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FIG. 5. d x 2_ y 2 (triangles) and s-wave (circles) pair-field correlations versus l/L. Those sin 
lations were carried out at $ = (see Eq. @). 
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FIG. 6. (a) n(k) at W/t = 0.6, U/t = 4 and (n) = 1. Lattices form L = 8 to L = 16 were 
considered, (b) same as (a) but for W/t = 0. The calculations in this figure were carried out at 
$ = (see Eq. (|)). 
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AE($, T) U/t = 4, W/t = 0.35, (n) = 1 




T/t T/t 



FIG. 7. (a) AE(&,T) as defined in equation (|3l]). (b) Vertex contribution of the equal time 
pair- field correlations as denned in equation (|36|). Here we use periodic boundary conditions in 
both lattice directions (i.e. $ = in Eq. @). 
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FIG. 8. One-electron density of states as a function of temperature and lattice size in the 
superconducting state at W/t = 0.35 and U/t = 4. L = 6: Figs, (a) — (c). L = 8: Figs, (d) — (/). 
L = 10: Figs. (5) - (m) 
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FIG. 9. (a) S(L/2, L/2) versus l/L for several values of W/t. The solid lines correspond to 
least square fits of the QMC data to the form l/L. Inset: S(L/2,L/2) versus l/L at W/t = 0.6. 
The solid is a least square fit to the form L~ a . (b) Staggered moment as obtained from (a) versus 



W/t. The data point at W/t = is taken from reference IE]. At Wtjt = 0.3, we were unable to 
distinguish m from zero within our statistical uncertainty. The solid line is a guide to the eye. The 
calculations in this figure were carried out at $ = (see Eq. (§)). 
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W/t = 0, (n) = 1 
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FIG. 10. (a) Uniform (b) staggered static spin susceptibility for the Hubbard model (i.e. 
W/t = 0) at U/t = 4. The solid line with no symbols in (a) corresponds to x(q = 0, = 0) at 
U = W = 0. The calculations were carried out at $ = (see Eq. (0)). 
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W/t = 0.35, U/t = 4, (n) = 1 
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FIG. 11. (a) Uniform (b) staggered static spin susceptibility in the superconducting state at 
W/t = 0.35 and U/t = 4. The calculations were carried out at $ = (see Eq. (||)). 
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FIG. 12. 5(Q, in the superconducting state at W/t = 0.35, U/t = 4 as a function of system 
size and temperature. For comparison, at the lowest considered temperature, T = O.lt, we have 
included the result for the Hubbard model (i.e. W/t = ) at U/t = 4. Here, $ = (see Eq. @). 
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FIG. 13. (a) 1/Ti versus T/Tkt in the superconducting state at W/t = 0.35. The magnetic 
scale is given by J ~ 0.5t. In the thermodynamic limit, Tkt ~ 0.2t and J is to first approximation 
size independent. Thus J/Tkt ~ 2.5 at L — > oo. (b). 1/Ti versus T for the antiferromagnetic 
Mott insulating state. Here, J ~ 0.25i. The calculations were carried out at $ = (see Eq. @). 
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FIG. 14. (a) Schematic phase diagram of the model H = Hjj + Hw — H*N at zero temperature 
and in the W — [i (see Eq. (|5|)). Here, ^ corresponds to the chemical potential and N is the particle 
number operator. Half-filling corresponds to fi = 0. At U/t = 4 , [i = 0, we estimate W c /t ~ 0.3 
and n c /t = 0.67 ± 0.015 Q. (b) Schematic phase diagram of H = Hjj + H\y in the T — W plane 
where T denotes the temperature. At U/t = 4 and W/t = 0.35 we estimate Tkt ~ 0.2i. 
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